Linear-response theory and lattice dynamics: a mufRn-tin-orbital approach. 
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A detailed description of a method for calculating static linear-response functions in the problem 
of lattice dynamics is presented. The method is based on density functional theory and it uses 
linear muffin-tin orbitals as a basis for representing first-order corrections to the one-electron wave 
functions. This makes it possible to greatly facilitate the treatment of the materials with localized 
orbitals. We derive variationally accurate expressions for the dynamical matrix. We also show 
that large incomplete-basis-set corrections to the first-order changes in the wave functions exist 
and can be explicitly calculated. Some usefuU hints on the k-space integration for metals and the 
self-consistency problem at long wavelengths are also given. As a test application we calculate 
bulk phonon dispersions in Si and find good agreement between our results and experiments. As 
another application, we calculate lattice dynamics of the transition-metal carbide NbC. The theory 
reproduces the major anomalies found experimentally in its phonon dispersions. The theory also 
predicts an anomalous behavior of the lowest transverse acoustic mode along the (CCO) direction. 
Most of the calculated frequencies agree within a few percent with those measured. 
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I. INTRODUCTION. 

Response of electrons to a static external field is one 
of the important characteristics of a solid which can 
be uniquely determined within density functional the- 
ory (D£'|T)|J. The use of the local density approximation 
(LDA)cHj in the linear-response problem has by now be- 
come a common and well established method for deter- 
mining various ground-state properties of real materials. 
These, first of all, include static dielectric, structural, and 
vibrational properties of semiconductors and insulators 
such as screening response to point charges and electric 
fields, effective charges and dielectric constants, as well 
as whole phonon spectraETij. Ab initio calculations of the 
wave-vector dependent lattice dynamical properti(^ for 
metallic systems have also recently been performedEiEJ. 
Among them, transition metals, their alloys and com- 
pound provide one of the fascinating areas in the study of 
phonons in crystal lattice. This is because, in addition to 
the richness and variety of structure of their phonon dis- 
persion curves, these materials also often exhibit lattice 
instabilities and relatively high-temperature (8K-23K) 
superconductivity, in which phonons play a fundamental 
role. Here, density functional based linear-response cal- 
culations provide an important first step in studying such 
phenomena as j^ectron-phonon interactions and trans- 
port propertiestHI which describe phonon-limited electri- 
cal and thermal resistivities, renormalization of specific 
heat (electron-mass enhancement) as well as supercon- 
ducting transition temperatures. These properties are 
connected with the real low-energy excitation spectrum 
of a metal. 

Initially, two methods have been developed to deal 
with the perturbations which break the periodicity of 
the original lattice. The first one, known as direct or su- 



percell approachll3, considers the perturbation with wave 
vector q which is periodical in the new supercell struc- 
ture. This is possible if the wave vector is commen- 
surate with the reciprocal lattice of the supercell and 
only tractable computationally if the size of the super- 
cell is not large. This limits the applications to high- 
symmetry wave vectors. The same technique can be 
applied to caJculate the interplanar force constants in 
direct spacell3. The dynamical matrix is found for any 
q using the Fourier transform provided the interatomic 
interactions of a solid are of short range. Despite of 
the severe computational-cost restrictions, the supercell 
approach has two important advantages: (i) the elec- 
tronic response and lattice dynamics can be studied using 
programs for self-consistent band-structure calculations 
which are standardly used in condense-matter physics, 
and, as a consequence, (ii), all non-linear-response coef- 
ficients are easily obtained. Note, however, that third- 
order non-linear coefficients can also be accessed within 
the linear-response approacht3 just like forces are found 
within the density-functional total-energy method. We 
shall discuss this statement in more detail later in this 
paper. 

The second method to deal with the perturbations 
is known as the perturbative approach. If the exter- 
nal perturbation is weak one can use standard perturba- 
tion theory and expand the first-order corrections to the 
one-electron wave functions in the unperturbed Bloch 
states of the original crystal. Usually it is done by in- 
troducing the so-called independent -particle polarizabil- 
ity function in terms of which the screened perturbation 
is found, by inverting the static dielectric matrix of a 
crystallla . Previously, due to a rapid progress made in 
the microscopic theory of the phonon spectra in free- 
electron-like metals through the development and ap- 



1 



plication of the pseudopotential technique, plane-wave 
representation was used for all the relevant quantities in 
the calculation. However, already in the case of covalent 
semiconductors with sufficiently weak pseudopotential, 
the convergency of the polarizability with respect to a 
number of plane waves becomes slow and there are only 
a few attepipts to compute phonon spectra within this 
frameworktj. The situation becomes worse for materials 
with localized orbitals. The most time-consuming step 
in this approach is connected with the problem of sum- 
mation over high-energy states which at least requires 
their calculation. Another problem is connected with the 
inversion of a large dielectric matrix. 

The above mentioned drawbacks of the perturba- 
tive approach have . beei n circumvented using the soliir 
state generalizationOEd of the Sternheimer methodE3. 
In this reformulated linear-response method, the first- 
order corrections to the unperturbed wave functions are 
found by solving the Sternheimer equation (which is 
the Schrodinger equation to linear order) directly with- 
out using the expansions over unperturbed states. This 
avoids the summation problem of the perturbation the- 
ory. The screening of the external perturbation is cal- 
culated self-consistently within DFT in close analogy 
with what is done in standard band-structure calcula- 
tions. This avoids the inversion problem. The present 
formulation is thus very efficient computationally which 
is demonstrated by an increasing number of applications 
to the problem of lattice dynamics in recent yearsHliJ. 

In order to solve the Sternheimer equation, one has to 
construct a rapidly-convergent basis set for representing 
the first-order perturbations. This is important because 
these corrections as well as the unperturbed wave func- 
tions oscillate in the core region. In the free-electron- 
like metals, broad-band semiconductors and insulators, 
this problem can be eliminated by the pseudopotential 
approximation and the majority of the a«glieations per- 
formed so far use plane-wave basis setsBflM. Unfortu- 
nately, with decreasing bandwidth, the plane-wave ex- 
pansion of the pseudo wave functions converges more 
slowly and it becomes less advantageous tojetnploy pseu- 
dopotentials. Indeed, until most recentlyailJ, the litera- 
ture contains no ab initio calculations of phonon disper- 
sions for transition-metal systems. 

In the present paper we describe an efficient all- 
electron generalization of the linear-response approach 
introduced in Refs. ^21 (A. brief report of this work 
has been published alreadyt£l). The first-order correc- 
tions are represented in terms of the muffin-tin (MT) 
basis sets which greatly facilitate the treatment of lo- 
calized valence wave functions. While the approach de- 
veloped in the paper is general and can be applied to 
any kind of localized-orbital representaticp. we use the 
linear-muffin-tin-orbital (LMTO) methodtj as a frame- 
work of such all-electron formulation. 

There are two problems addressed in this paper which 
are connected with the use of MT basis functions in 
a linear-response method. The first problem concerns 



the construction of a variational solution of the Stern- 
heimer equation. This is necessary because the unper- 
turbed energy bands and wave functions are obtained 
within the LMTO method by applying the Rayleigh- 
Ritz variational principle. They are, therefore, not ex- 



act solutions of the one-pS 
As first shown by PulayE 



Jectron Schrodinger equation. 
3, the use of variational solu- 
tions gives rise to the incomplete-basis-set (IBS) cor- 
rections in force calculations. The IBS corrections must 
be carefully accaunted for to get accurate forces in the 
LMTO methodnJ. These corrections also exist and must 
be taken into account when calculating the first-order 
changes in the wave functions and the dynamical matrix 
within the linear-response approach. The other problem 
is connected with the change in the basis functions due to 
the perturbation. Since the one-electron wave function 
in the LMTO method is represented by the expansion co- 
efficients in the MT basis set, under static external per- 
turbation, such as the displacement of a nucleus from its 
equilibrium position, the change in the wave function will 
be described by both the change in the expansion coeffi- 
cients and the change in the basis set. The contribution 
from the change in the basis set is important because the 
original basis set is tailored to the one-electron poten- 
tial and must therefore be reconstructed to account for 
the specifics of the perturbation. It should be noted that 
this contribution is not taken into account in the stan- 
dard perturbation theory, where only the change in the 
expansion coefficients is taken into account. 

The rest of the paper is organized as follows. The 
variational formulation of the linear-response approach 
is described in Section II. Implementation of the theory 
in the framework of the LMTO method is described in 
Section III. Application of the method to phonon spectra 
in Si and NbC is given in Section IV. Section V concludes 
the paper. 



II. THEORY. 

a. Density-functional linear response. 

Within density functional theory, the problem of cal- 
culating the lattice dynamics essentially amounts to find- 
ing the change in the electronic charge density induced 
by the presence of a phonon with wave vector q. Con- 
sider a lattice with a few atoms in the unit cell given by 
the positions R -I- t, where R are the basis vectors and 
t are the translation vectors. Suppose that the atoms 
are displaced from their equilibrium positions by a small 
amount: 



6tii = SARexp{+iqt) + SA*j^exp{—iqt), 



(1) 



where SAf; is a complex polarization vector and q is the 
phonon wave vector. The presence of such displacement 
field changes the bare Coulomb potential as follows 
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R.t 

where Zfj are the nuclei charges. By expanding this ex- 
pression to first order in the displacements we obtain that 
the crystal is perturbed by the static external field: 

R t ' I 

y SA*j,y e-^^'W , y , , (3) 
^ r-R-t ' ^ ' 

R t ' I 

which is represented as a superposition of two travelling 
waves with wave vectors +q and — q, i.e. 

SVext{r)^}_^SAR +}_^SAj, . (4) 

R R 

To shorten the notations, we will omit (5R from this 
definition and, therefore, SVext = J2^-^^^^ext + 
^ (SA*i5~T4a;t- Both components S^Vext and 5~Vext 
are hermitian: [(J^Ve^jt]* = S^V^xt and translate like 
Bloch waves in the original crystal: S^Vexti^ + t) = 
exp{±iqt)S^Vext{ic) ■ 

The first-order change in the charge density, 6p, in- 
duced by the perturbation (H) is represented in the same 
form as (||), i.e. 6p = Y^-^^^^ P + J2 ^ ^~ P and it 
is expressed in terms of the one-electron wave functions 
ipu_j and their first-order corrections S'^ipi^j and 5~^-kj as 
follows 

5^P = J2 hjiS^^wji^^j + KjS^i'wj), (5) 

where /kj are the occupation numbers, k lies in the first 
Brillouin zone, and j is the band index. The first-order 
correction S^ipkj = |(5^kj) = (S^ipkj)* is a Bloch wave 
with vector k±q and it is the solution of the Sternheimer 
equation, which is the Schrodinger equation to linear or- 
der: 

+ Veff - ey,,)\6^kj) + 6^V,fj\kj) - 0. (6) 

Here, 14// is the effective DFT potential and S^Veff is 
the change in the potential which is the external pertur- 
bation screened by the induced charge density: 

^-ye// = ^-K. + e^/^ + ^^-p, (7) 

where the exchange-correlation effects are taken into ac- 
count in the local density approximation. In Eq. (|^) we 
have dropped the term with the first-order corrections 
to the one-electron energies: (5^ekj = {kj\S'^Vef f\k.j) , 
which are equal to zero if q ^ 0. The Eqs.(||-(|^) must be 
solved self-consistently, i.e. (i) one has to solve (^) with 
the external perturbation S^Vext, or the one screened by 



some guessed S^p, (ii) find the induced charge density ac- 
cording to (H) and, (iii), find new S^Veff after (0). Steps 
(i)-(iii) are repeated until input and output 6^p coincide 
within a required accuracy. This is much analogously to 
finding the unperturbed quantities ipi^j , p and 14/ / in 
standard band-structure calculations. 

b. Dynamical matrix. 

We must now solve two problems in order to calculate 
the lattice dynamics: to develop a method for solving the 
Sternheimer equation (^, and to find an expression for 
the dynamical matrix. The general strategy employed 
in the following is to consider the dynamical matrix as 
a functional of the first-order perturbations. Expanding 
S^4'k.j in terms of the MT-basis functions will lead, under 
the stationarity condition, to a set of matrix equations 
which represent a variational solution to Eq. (H) . 

The variational formulation of the linear-response 
problem is required because the original states i/'kj are 
normally found not as exact but variational solutions of 
the one-electron Schrodinger equation. In an all-electron 
method such as the LMTO method, the wave function is 
expanded in terms of the basis set |Xq): 

i^^j=J2\Xa)A^J- (8) 

a 

The total energy is then considered as a functional of 
only the expansion coefficients , which are found by 
applying the Rayleigh-Ritz variational principle. This 
leads to the following matrix eigenvalue problem: 

E(^^l - + ^^// - ek.lx^)^^^' = 0, (9) 

a 

which, in particular, allows all non-spherical terms in the 
potential to be taken explicitly into account . 

In the problem of lattice dynamics the second-order 
change in the total energy must be found. It is obtained 
by expanding the total energy with respect to the change 
in the external potential (nuclei displacements) up to sec- 
ond order, i.e. E = Eq + 5^^^ E + 5^'^^E. The first-order 
change S^^^'E vanishes if the lattice is in the equilibrium 
and the second-order change is expressed via the dynam- 
ical matrix of a solid: 

S^'^E^^ ^ A^4(q)M«v"5A^,^ + c.c., (10) 

R' ^j,' Rfi 

where we assumed that the nuclear displacements have 
the form (|l]) and where {p} denote directions of the po- 
larization. A formula for dynamical matrix A is obtained 
by varying a density-functional expression for the total 
energy. It is given by 
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A^4(q)^A = 

J2 /k,('5+'5-kj + S-S+kj\ - + Veff - 6k,|kj) 
J2 /k,2(5+kj| - V2 + - ek,|5+kj) + 



kj 



(11) 



(We omit indexes R'n'Rp for simplicity.) Here j^^J^kj) 
denote second-order changes of the wave functions, and 
d~^6~Vext is the second-order change in the bare Coulomb 
potential due to the displacements of nuclei. The latter 
is given by 



0^0 Vext = Off/R > V^'V^-j — 

|r — R - 



(12) 



Note that, when expanding Eq. (g), it is only necessary 
to keep the terms periodical in the original lattice, be- 
cause only these give rise to a non-zero contribution to 
the integral with p. Consequently, all the contributions 
with the phase factor ea;p(±2zqt) can be neglected. The 
same is true for the second-order terms: IS^S^kj) are 
the functions with wave vector k since they contribute 
to the matrix elements with |kj) only. (Therefore, the 
operator can be considered as a variation of a Bloch 
wave, whereby ±q gets added to its wave vector.) 

The first term in expression (|l^) is not zero if the un- 
perturbed states are approximate solutions found from 
the eigenvalue problem (g. If, on the other hand, one 
neglects it, and performs a variation of A with respect 
to the first-order corrections, the self-consistent linear- 
response equations (^), (|^) will be recovered. The ex- 
pression (^) is variational with respect to the first-order 
changes in the wave functions just like the unperturbed 
total energy is variational with respect to the unper- 
turbed states |kj). This property of the density func- 
tional follows from the Hohenberg-Kohn-Sham varia- 
tional principle. Eq. ( |ll|) is directly interpreted as the 
electronic contribution to the dynamical matrix. (The 
part connected with change in the bare Coulomb energy 
of the nuclei is evaluated trivially by the Ewald tech- 
nique.) Due to stationarity of this functional the calcu- 
lation of the dynamical matrix is very accurate: while 
the first-order changes in the wave functions and the 
charge densities are only variationally accurate, the er- 
ror will be of second order with respect to the error in 
|(5^kj). In particular, the convergence of the dynamical 
matrix during the iterations towards self-consistency is 
much faster than the convergence of the induced charge 
density. At its minimum, expression (11) contains no 
second and third terms and it may therefore be inter- 
preted as the Hellmann-Feyman result (last two terms) 
plus incomplete-basis-set correction [the first contribu- 



tion]. The latter is the analog of the "Pulay force" known 
from force calculations. 

The discussed variational properties of the dynamical 
matrix are indeed not unique and represent a particu- 
lar case of the powerful "2n + 1" theorem of perturba- 
tion theory: the knowledge of the perturbations in the 
wave functions up to (n)th order allows one to find the 
(2n -I- l)th correction to the eigenenergy. A generaliza- 
tion to arbitrary order of perturbation within the density 
functional theory as well as variational properties of even 
derivatives of the total energy were discussed in Ref. |l^. 
A direct minimization of the dynamical matrix by the 
conjugaie-gradient method has also been demonstrated 
recentlyQ. It worth to point out that the knowledge of 
the first-order corrections to the wave functions allows us 
to consider changes in the total-energy up to third order, 
in the same way as the zeroth order unperturbed states 
allow calculating such first-order derivatives as, for in- 
stance, forces. Consequently, third-order unharmonicity 
constants, Gruneizen parameters and other non-linear 
coefficients are, in principle, easily accessed within the 
linear-response formalism. 

c. First-order corrections. 



We now turn out to the construction of the basis func- 
tions which represent the first-order perturbations. Since 
the unperturbed state is given by expansion (g) , the first- 
order change |(5^kj') must include both the change S^A^^ 
in the expansion coefhcients as well as the change |(5*xJ^) 
in the MT basis set, i.e. 

\S^kj) = J2{\xl+'^)S^Al= + IS^x^JAl^}- (13) 



Since |(5^kj) is a Bloch state with wave vector k ± q so 
are IXq^'') a-nd If^^Xa) ■ The first function is the original 
linear MT orbital of wave vector k±q and the second one 
is the change in the MT orbital due to the movements of 
atoms. In Section III we will give detailed formulae for 
the change in the basis functions. Here we note that since 
the original basis |Xq) is a Bloch sum of atom-centered lo- 
calized orbitals, the important contribution to the change 
in the Bloch sum is connected with the rigid movement 
of these orbitals due to the rigid movement of the poten- 
tials for displaced atoms. The expansion (|3|) is rapidly 
convergent because the basis can be tailored to 

the perturbation just like the basis |Xq) is tailored to the 
unperturbed potential. Eq. (^3|) can be interpreted as an 
expansion of |(5^kj) in terms of Ix^"'''*) in the local coor- 
dinate system displaced with the atom; the convergence 
with respect to the number of orbitals per atom must be 
about the same as for the unperturbed state. This is in 
contrast with the expression of the standard perturba- 
tion theory where, for the expansion of \S^hj), only the 
change in the coefficients A^^ is taken into account [first 
contribution to (Ohl. 
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We shall now consider the second-order changes 
I^^J^kj) which appear in expression ([ll| ) for the 
incomplete-basis-set corrections to the dynamical ma- 
trix. By performing the variation of the expansion (||) to 
second order we obtain: 

+ \S^xl'''')S^A^J + \S^S^xl)Al^, (14) 

where S^S^A^^ and |<^*<^^Xq) are the second-order 
changes in the expansion coefficients and the basis func- 
tions respectively. By inserting ( p^ in the first term of 
( |rT| ) one sees that the second-order changes S^S^A^^ do 
not contribute because they enter as coefficients to the 
unperturbed basis functions and 



S^S^A^^X^a\H-eu,\kj)^Q. 



(15) 



(Here, H = -V^ + Veff.) 

The absence of the coefficients 6^6^ A^^ in our for- 
mulation of the problem has an important consequence: 
since |(5^(5^kj) has only the unknown contribution from 
the first-order changes in A^^ and since the Hilbert space 
{|x),|^x)} of the basis functions is fixed, we see that 
the variational freedom of the functional ([ll| ) is provided 
only by the coefficients S^A^K This is again in close 
analogy to that in band-structure calculations the vari- 
ational freedom of the total energy is provided only by 
the unperturbed coefficients A^^ . In the total-energy cal- 
culations this has the consequence when calculating the 
forces: due to the stationarity condition the force formula 
does not contains any first-order derivatives in A^^ . In 
the dynamical-matrix calculation this will have the same 
consequence when calculating third-order non-linear co- 
efficients: the corresponding formulae will not contain 
any second- and third-order derivatives of ^J^-' and, thus, 
can be explicitly evaluated from only the knowledge of 
S^A^^ . Note however, that together with the matrix el- 
ements containing |x) Ji^x) J'^^^^x); a contribution from 
third-order changes in the basis sets must be taken into 
account. 

We shall now derive the equations for the first-order 
changes in the expansion coefficients. This is done by 
minimization of A with respect to 6^A^^ . We obtain: 

a a 

(16) 

This linear system of equations is, in fact, a variation of 
the original eigenvalue problem (|9|). It determines the 
position of the minimum of A in the space of the coeffi- 
cients S^A^^ , and non of the second-order changes, such 
as |i5='=i5''';^;^), aflFect it. The functions \S'^S^Xa)y on the 
other hand, define the value of A itself in its minimum 



and must be taken into account in the evaluation of the 
dynamical matrix. 

We must now solve equation (p^). This equation in- 
volves only the occupied states of the unperturbed sys- 
tem, which are necessary for constructing the induced 
charge density according to (||). It may be solved us- 
ing an iterative algorithm with the number of opera- 
tions proportional to Nband x ^basis^ where N^and is a 
number of filled bands and Nbasis is a number of the 
basis functions used for representing the unperturbed 
wave functions and their first-order corrections. This 
scheme is advantageous when using the LAPW or plane- 
wave pseudopotential methods where the conventional 
matrix-diagonalization algorithms represent the most 
time-consuming step which scales as the cube of the size 
of the basis. The LMTO method, on the other hand, has 
a small basis and the inversion of the matrix {x^'^\H — 
EkjlXa^'') required for solving Eq. ([l^), can easily be 
performed because its eigenvalues are eiciqj' — Ckj and 
eigenvectors are A}^'^ , where f = l,Nbasis- Because 
of the minimal size of the basis in the LMTO method, 
it is not a time-consuming step to find eigenvalues e-^j 
and eigenvectors A^^ for all energy bands {— Nbasis) at 
some grid of wave vectors k. This is independent of the 
phonon mode and therefore needs to be done only once. 
We therefore use the original eigenstates for the matrix 
inversion. The result for 6^A^^ is then substituted into 
( p^ ) which gives the final expression for |(5^kj) in the 
form: 

i<5±k,> y: i^^^^)^!^^' + E ^^^^ X 

{{k±qj'\H ~ eu,\J2S^xlAl^) + 

a 

{Y,S^xl''''Al^'"'\H-e^,\kj) + 

a 

{k±qj'\S^Veff\kj)}. (17) 

This formula has a simple physical meaning. The first 
three terms containing \Sx) appear because of the use 
of variational solutions. They can be interpreted as 
incomplete-basis-set corrections to the last term (the 
one with S^Veff), which has the form of the standard 
perturbation theory. If all unperturbed states are exact 
and they represent mathematically a complete basis set, 
then the first and second terms in ( [l7| ) cancel, the third 
term vanishes and the standard pcrturbative formula is 
recovered. However, if this is not the case, the use of the 
functions \Sx) in the basis greatly reduces the number of 
states |k ± qj') needed to reach the convergence in (|T^). 
Namely, following the above derivation, the summation 
in the last three terms is over Nbasis energy states, i.e. 
over the size of the basis for the unperturbed system. 
To illustrate the advantage of this formula we consider 
the so-called acoustic sum rule ( ASR) : suppose all atoms 
are displaced in the same direction by a small amount. 
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The change in charge density induced by the rigid move- 
ment We// of the potential will be equal to Vp. Within 
standard perturbation theory one obtains: 



\Skj) = ^ \kj 



.,0^f\VVeff\kj) 



^|k/)(k/|V|kj)=V|kj). 



(18) 



The latter equality can only be obtained if the states 
|kj') represent a mathematically complete basis set. This 
is not case in the LMTO method which employs a min- 
imal basis set to reproduce the energy bands and wave 
functions within a certain energy window. On the other 
hand, within the minimal basis set the ASR can be triv- 
ially satisfied if one uses the expression for the first- 



order corrections: here, by construction \Sxo 



while the last three contributions vanish. (This is so be- 
cause they are combined into the integral from a gradient 
of the periodic function: Vjf/'kj/ {H — ekjOV'kj}, which is, 
by definition, equal to zero.) 



The unoccupied states in the expression (17) should 



not be considered as real excitation energies and wave 
functions. Let us consider the induced charge density as 
a ground-state property of both perturbed and unper- 
turbed systems. In both cases only the occupied states 
must be well reproduced, the excited states can, in prin- 
ciple, be arbitrary. The LMTO and LAPW methods are 
very suitable for this purpose: they are fast and accu- 
rate within a certain energy window, which is achieved 
by expanding the basis functions of the original KKR 
and augmented-plane-wave (APW) methods by Taylor 
series around some energies at the centers of inter- 



est. The states |k ± qj' 



|17|) are the eigenstates 



of the Hamiltonian matrix {Xp^'^\H\'X^'^) which is it- 
self constructed to reproduce the occupied energy bands 
well. This is the energy window of interest and all cen- 
ters of linearization are in this window. In the KKR 
and APW methods the states |k ± qj') have the follow- 
ing meaning: since the KKR (APW) energy bands and 
eigenvectors are the eigenstates of the LMTO (LAPW) 
Hamiltonian {')^{^u)\H\}^{e^)) with = Ckj, the states 
|k ± qj') in ( ll7| ) will be the eigenstates of the Hamil- 
tonian {Xp^'^{^\ij)\H\^^^{e]ij)) and only those bands 
Ekiqj' with energy near ekj will be described correctly. 
In this case, finding |(5^kj) requires the knowledge of 
this auxiliary spectrum for every occupied energy ekj. 
We thus finally conclude that the excited states are not 
to be interpreted as the exact ones, only the knowledge of 
occupied energy hands is necessary in our linear-response 
formulation. 



III. IMPLEMENTATION. 

In this section, an extension of the linear muflin-tin 
orbital method for linear-response calculations is de- 



scribed. We shall first review the full-potential LMTO 
method, which is used as the framework in this im- 
plementation. Then, the problem of constructing the 
changes in the MT-orbitals due to the atomic move- 
ments is discussed. Other problems considered are the 
Brillouin-zone integration for metallic systems and the 
self-consistency at long wavelengths where the Coulomb 
singularity Air/q'^ makes the standard mixing schemes 
computationally inefhcient. 

a. Full potential LMTO method. 

We first review the LMTO method, which solves the 
original Schrodinger equation. The space is partitioned 
into the non overlapping (or slightly overlapping) muffin- 
tin spheres sr surrounding every atom and the remain- 
ing interstitial region f2i„t. Within the spheres, the basis 
functions are represented in terms of numerical solutions 
of the radial Schrodinger equation for the spherical part 
of the potential multiplied by spherical harmonics as well 
as their energy derivatives taken at some set of energies 
ei, at the centers of interest. In the interstitial region, 
where the potential is essentially flat, the basis functions 
are spherical waves taken as the solutions of Helmholtz's 
equation: (— — e)/(r, e) = with some fixed value 
of the average kinetic energy e = k^. In particular, in 
the standard LMTO jaethod using the atomic-sphere 
approximation (ASA)Ej, the approximation = is 
chosen. In the extensions of the LMTO method for a 
potential of arbitrary shape (full potential), a multiple- 
kappa basis setcj is normally used in order to increase 
the variational freedom of the basis functiona_.while re- 
cent developments of a new LMTO techniquecJ promise 
to avoid this problem. 

The general strategy for including the full-potential 
terms in the calculation is the use of the variational 
principle. A few different techniques have been devel- 
oped for taking the non-spherical corrections into ac- 
count in the framework of the LMTO method. They 
include Foucicj; transforms of the LMTOs in the intersti- 
tial regionEaEj, one-ceiiifir spherical-harmonics expan- 
sions within atoraic cell£3, interpolations in terms of the 
Hankel functionsEJ as well as direct calculations of. the 
charge density in the . tijght-binding representationEiJ. In 
two of these schemesOE3 the treatment of open struc- 
tures such as, e.g. the diamond structure is complicated 
and interstitial spheres are usually placed between the 
atomic spheres. In the dynamical-matrix calculation it 
is inconvenient to use interstitial spheres because they 
lead to artificial degrees of freedom for the lattice dy- 
namics. Therefore we will develop the linear-response 
LMTO technique using the plane-wave Fourier represen- 
tation. This allows us to apply the method for such ma- 
terials as Si and NbC without interstitial sphexps..|J^ote, 
however, that in our previous applicational3il3'E3 for 
BCC and PCG metals, atomic-cell spherical-harmonic 
expansionsEa were used. 
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Consider the so-called envelope function, which is a 
singular Hankel function, 

K^RLirR - t) = K^RiilvR - t\)i^Yim{rR - t), (19) 

centered at site R + t and with an energy e = nf,. Yim de- 
notes a complex spherical harmonic with the phase con- 
vention after Ref. Inside any other site R' -I- t' the 
Hankel function can be represented as an expansion in 
terms of the Bessel functions, JnR'L'i^R' — t'), i.e 



K^rlIyr 



L' 



t) 

JkR'L' 



{vw - t'hRn'SR'L'RLit' - t, k), (20) 



where jri = 1/sr{21 + 1) and Sr'l'rl{^, k) are the struc- 
ture constants in real space. Note that, while the index 
L enumerating the basis functions usually runs only over 
s,p, and d states, the sum over L' in this expression must 
include higher angular momenta. Normally goes up to 
6 — 8. For convenience, we use the following prefactors 
in the definitions of the spherical functions: 



K,,Ri{\rR\) = 



JkRiUvrI 



{'iSRy+^ 



(2Z + 1)!! 



(tSfl)' 



(21) 



(22) 



where hi = ji — ini , ji , and ni are the spherical Han- 
kel, Bessel, and Neumann functions, respectively. The 
expression for the structure constants is then 



5™(t,.) = (^)''"V^) 

\ w / \ w / 



E 



-4Trw{2l" - 1)!! 
■;f (2/'- 1)!!(2Z- 1)!! 



i+i'-i" 



7^,^P-(|t-R' + R|)(-z)' y£„(t-R' + R), (23) 

where w is the average Wigner-Seitz radius and the Han- 
kel function K^^^i is defined with w instead of sr in ex- 
pression (|2l]). The Gaunt coefficients C^j^i are defined 
by the integral: 



LL' 



(24) 



We now consider a Bloch sum of the Hankel functions 
(|l9|), centered at different sites, which, inside the MT- 
sphere at R', is represented by the expansions in the 
Bessel functions: 



ikt 



kRL 



(l-fl - t) = 



K^,rl{vr)5r'r - ^ Jf,R,L'{rR'hR'i'S%i^,jiL{K), (25) 

L' 

where S^, ^ , j^^Xk) denotes the lattice sum of the struc- 
ture constants (|2^). The linear MT-orbitals \x^rl) 



are now obtained by augmenting the spherical functions 
KkRl,JkRL in all MT spheres by numerical radial hmc- 
tions ^^rl.Hrl- 



XkrA'^B') = 



K 



'RLirR)SR'R - E <f«fl'L'(ri?')7i?'i'5|,i,«i(K). (26) 
r I 



The functions ^^r^ 



kRL 



are the linear combinations 



of the solutions ipRLiJ^R, ^i^kRi) = 4>kRL to the radial 
Schrodinger equation as well as their energy derivatives 
4>RLirR,e^KRi) = 4>kRL taken at the energies c^kRI- In 
the interstitial region, the linear MT orbitals are repre- 
sented as multicenter expansions [left-hand side of Eq. 
(^)]. In order to calculate the interstitial-potential ma- 
trix elements and represent the charge density, we use 
the Fourier transform of the LMTOs in the interstitial 
region. It is impossible to consider the Fourier transform 
of the expression (|2^) directly because of the singulari- 
ties in the Hankel functions. On the other hand, since 
this representation will be used for the description of the 
basis functions only within Qint, we can substitute the 
divergent part of the Hankel function by a smooth func- 
tion for vr < Sr. This regular function is defined in the 
Appendix and it is denoted as Kj^rl. We thus introduce 
a pseudoLMTO IXkAl) defined in all space as follows: 



XkRL 



(r) = ^e*'^*if«x(r«-t) = 

t 

= ^XKKL(k-^G)e^(''+°)^ 



(27) 



G 



which is identical with the true sum ( p5| ) in the intersti- 
tial region. 

The charge density and the potential have a dual repre- 
sentation: spherical-harmonic expansions inside the MT- 
spheres and plane-wave expansions in the interstitial re- 
gion. This is usually done by introducing a smooth pseu- 
docharge density_p in all space defined in terms of the 
pseudoLMTOs (|27|). The pseudodensity coincides with 
the true density when r g ^int ■ In this way, the solu- 
tion of the Poisson equation is straightforward and can 
be done along the lines developed in Ref. 33 . In practical 
applications we have also used the technique described in 
the Appendix for the Fourier transform of the Coulomb 
interactions and for the construction of auxiliary densi- 
ties. The exchange-correlation potential is found using 
the fast Fourier transform and the interstitial-potential 
matrix elements are explicitly evaluated. 

b. Changes in the linear mufRn-tin orbitals. 

We shall now discuss the linear-response calculation. 
Small displacements of atoms from their equilibrium po- 
sitions defined by expression (||) lead to the change in the 
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Bloch sum of the atom-centered (pseudo) Hankel func- 
tions (^tI). Because of the exphcit dependence of the 
basis functions |Xkj?x) '^^ the atomic positions R, here 
and in the following the displaced atoms will be denoted 
by the index R. The change in the Bloch sum of the MT 
orbitals can be found analogously to the change in the 
external potential in Eqs.(^-(P. As a result, we con- 
sider two travelling waves with wave vectors k -|- q and 
k — q, i.e. 



5Ru 



-Srr Y.'ik±q + G)^ X x«flL (k ± q G)e'(''±'^+°)^ 
G 

(28) 

which represent the change of the basis functions in the 
interstitial region or the change of the pseudoLMTOs in 
the whole space. Here we have restored the original nota- 
tions: 5'^^rl{^)I5R^j, = ^'^ji^-RLi^)- We also introduce 
a spherical coordinate systemtHI: 



R 



(29) 



4'kRL is a result of both the rigid movement of the spher- 
ical component of the potential and the change in the 
shape of the spherical component. In the following, it 
is convenient to treat the rigid movements of the poten- 
tial within the MT-sphere centered at R separately, i.e. 
represent the total change in the form: 



SRa 



St)Veff{rR) 



6Ru 



(32) 



where the notation stands for the "soft" change, i.e. 
the variation connected with the change in the shape 
of the function. The functions 5^(j)K,RL/SR^ are repre- 
sented in a form similar to (^^, i.e. 5^4>KRL/5Rfj, = 
-Srr^ ii4>KRL + 5^s)'l>KRL/5Rfj,, where the last (soft) con- 
tribution is found by solving the radial Sternheimer equa- 
tion: 




which is connected to the Cartesian system as follows: 
R~i — +{Rx — iRy)/V2, Ro = Rz,R+i = —{Rx + 
iRy)/V2. The reason is that, in the spherical coordinates, 
the operation on a product of a radial function /(r) 
multiplied by the spherical harmonic takes the simple 
form: 



V^/(r)rz„,(r) - 



I 



fir) li+i™+M(r) + 



(30) 



We shall now find variation of the basis functions inside 
the MT spheres. In the sphere R', the original LMTO is 
defined in the expression (^). Its change must include 
both the changes in the numerical radial functions and 
the change in the structure constants: 



E 



5Ru 



5R,. 



R'R 



lR'vS\,L,jiL{K) 



5Ru 



(31) 



The change in the numerical functions contains two con- 
tributions. Since ^kRL constructed from the 
solutions of the radial Schrodinger equation and their en- 
ergy derivatives, 4>k,rl and tp^RL, the change in c/j^rl and 



The superscript " SPH" here denotes the spherical com- 
ponent of the potential and the perturbation. It is, in 
principle, not a problem to take all non-spherical terms 
of the perturbation into account. If this is done, the 
first-order changes in the radial functions are no longer 
given by a single spherical harmonic but as an expansion 
in Yim- One obtains an uncoupled system of radial equa- 
tions, which can easily be solvedo. However, in the prob- 
lem of lattice dynamics the change in (fi^RL and (puRL due 
to the change in the shape of the spherical component of 
the potential is small. This is because the motions of 
atoms mainly distort the dipole part of the potential. 
If the change in the shape of the spherical component 
can be described by some constant shift of the energy, it 
may be cancelled by appropriate choice of the change 
Sfg-jiuKRi/SRti in the energies ei^^Ri- This cancellation 
can, for instance, be obtained by finding S^^-^e^i^m/SR^, 
with fixed logarithmic derivatives D^t^m . (The deriva- 
tives D^i^m are evaluated at the occupied centers of grav- 
ities of the bands for the unperturbed crystal.) We thus 
see that the influence of the constant shifts to the change 
in the basis set can be eliminated and, therefore, one 
can neglect by the contribution S^^^fjif^nL/SR^ in practi- 
cal calculations. The accuracy of this approximation is 
quite good which has already been confirmed by good 
agreement between total energy j-and force calculations 
with the original LMTO methocuJ where the same ap- 
proximation was used for deriving the force formula. 

We now give the formula for the change in the struc- 
ture constants which enters Eq. (|3l|). It is expressed 
as the difference between the gradients of the structure 
constants for wave vectors k and k ± q, i.e 
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^Rr"^ ^i^a'L'RLi.'^)- 



(34) 



The gradient is with respect to R' — R. From (po|), using 
the recursion relations for the Hankel functions, it fol- 
lows that the change in the structure constants can be 
expressed in terms of the structure constants: 



^ fi^R'l'm'Rlmi'^) 



K^Sr> 



2 ^l' — Im' —fil'm' ^21' -^^ R'l' 



2 I' -\-lm' — fil' m' 



{21' + 1) 



SR' 



ck 

^R'l' + lm' 



lm'-/ifl/in('*) + 

{k). (35) 



Here, the left index of the structure constants has 
changed to /' ± 1, m' — but the right index Zm remains 
the same. An analogous formula exist in which the right 
index change is Z ± 1 , m + in and the left index is un- 
changed. 

It is seen that the change in the MT orbital (^ij) can 
be represented as a rigid part, a small soft part and a 
contribution from the change in the structure constants: 



SRu 



Sts)X^RLi^R' 

SRu 



E*««'^'('-^')7i.,'.^-%^. (36) 



It is convenient to separate the rigid part since it gives 
rise to a rigid contribution in the electronic response: 



SRu 



^ -SRR^f.p{rR) 



SRu 



(37) 



Since the induced charge densit y (p7|) has the same form 
as the change in the potential (|3^)7 we need not calcu- 
late the gradients of the charge density and the potential. 
This is important since these gradients are huge in the 
core region, which could result in large numerical errors. 
The second term in ( ^ is S'^^s^x^RLi^R')/SRti- It is con- 
structed from the changes S^^s^cj)njiL/SRi_i and their energy 
derivatives which are numerically small. This function is 
exactly equal to zero together with its first-order radial 
derivative at the sphere sri. It translates like a Bloch 
wave with vector k ± q because the original form of one- 
center expansion ( |26| ) translates with wave vector k while 
the first-order changes S'^^s^(j)i^jiL/ SRf^ translate, like the 
perturbation, with wave vector ±q. The whole expres- 
sion (36) also translates with wave vector k ± q and fits 
into the multicenter expansion of the change in the ba- 
sis set in the interstitial region [formula (|28|)]. However, 
since the original LMTOs are continuous and only dif- 
ferentiable to first order at the boundaries of the MT 
spheres, the matching of the change in the basis set is 



only continuous but not differentiable. This, in princi- 
ple, leads to a kink in the change of the charge density. 
However, it does not have any effect in the calculation 
of the dynamical matrix if the latter is compared with 
the second-order derivative of the total energy derived 
from the frozen-phonon supercell calculation. This is so 
because the extension of the LMTO method described 
here is just an analytical version of the finite-difference 
approach employed in the supercell technique. When ap- 
plied to the same problem, the results of both approaches 
have to be the same except for the errors introduced by 
taking finite differences. This concerns the comparison of 
not only the dynamical matrix and the phonon frequen- 
cies, but also the changes in the basis set, the expansion 
coefficients, the charge densities and in all other quanti- 
ties which can be obtained by the frozen-phonon LMTO 
technique. 

We now turn to the problem of calculating the change 
in the expansion coefficients ^^'^i, which are necessary 
to compute the first-order corrections according to (p^). 
From expression (|l^, the change S^^Jj^^/SR^ IS given 
by 



5^^RL 

SRu 



E 



/ikiqj' 

^kRL 

,5±/fk±qj'kj 

sE, 



jiQkicLj'kj 



(38) 



and it is expressed in terms of the change in the hamil- 
tonian and the overlap matrices. Here, the change in the 
matrix elements of the hamiltonian is given by the band 
representation: 



sWu 



/^\^±cij' * S^ H K,' R' L' kRL ^kj 
2^ ^K'R/L' 



k' R' L' K.RL 



SRu 



kRL 



(39) 



and a similar formula holds for the matrix elements of the 
overlap integral. In the original, {kRL}, representation 
the changes S^H^iRiL'^Rh/SRfj, aiid S^ O r> l' kRL / SRn 
are readily computed using the formulae ( |28| ) and ( |3q ) 
for the first-order changes in the basis set. It is indeed 
even more advantageous to find the corresponding formu- 
lae by directly varying the expressions for the hamilto- 
nian and the overlap matrices. This avoids the problem of 
combining the contributions with the gradients of numer- 
ical radial functions to the surface integrals. One point 
about calculating the change in the interstitial kinetic- 
energy matrix elements and the interstitial overlap in- 
tegrals is worth noticing. Since these matrix elements 
contain energy-derivative of the structure constants, the 
change in these matrix elements will contain the change 
in this energy derivative. The corresponding formula can 
be found by taking the derivative with respect to in 
the expressions (Q) and (|3^). 

Another problem is to find second-order changes in the 
LMTO basis functions as well as second-order variations 
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in the Hamiltonian and the overlap matrices. They are 
necessary for computing the incomplete-basis-set correc- 
tions in (|ll|) for the dynamical matrix. In the interstitial 
region the second-order change in the pseudoLMTOs is 
simply given by: 



G 



(40) 



Inside the MT spheres the expression is more compli- 
cated, but can be found straightforwardly by perform- 
ing one more variation /SRfj,' of expression (^TJ) for 
the first-order change. It will contain second-order 
changes in the numerical radial functions and second- 
order changes in the structure constants as well as differ- 
ent products of the first-order changes in these quanti- 
ties. The second-order changes in the structure constants 
are given by 



R'L'RL 



in) 



s-s 



k-l-q 

R'L'RL 



sK' 



(41) 



This is obtained from the expression ( p^ ) and is 
expressed via the difference between the gradients of the 
structure constants for the wave vectors k and k — q, 
while (5~ S^^"^ is the difference between VS for the wave 
vectors k-l-q and k. Alternatively, the expression ( ^ ) 
can be found by first considering the expression for the 
structure constants in the supercell and then, assum- 
ing the form ([^) for the atomic displacements, transfer- 
ring the supercell expression to the original structure. 
The second-order gradients V^' S are calculated us- 



ing (|3^) and they are again the structure constants with 
the left index changed to /' ± 2, m' — fi' — fj, and unchanged 
right index. Analogously, they can be expressed in terms 
of the structure constants of the same left index I'm' and 
the right index: I ±2,m + fi' + fj,. 

The second-order changes in the numerical radial func- 
tions must also be calculated. They contain contributions 
V^'V^^f^"^ due to the rigid movement of the spherical 
part of the potential to second order, changes due to the 
rigid movements of the first-order variations in the shape 
of the spherical part (rigid movement of the soft part), 
as well as the contributions arising from the change in 
the shape of the spherical part to second order (second- 
order soft part). As we discussed above, one can neglect 
by the infiuence of the change in the shape of Vfj:^^ to 
the change in the basis. Therefore, we must only keep 
the rigid contributions described by Vp'V^$^^. 



c. Brillouin-zone integrals. 



After computing the first-order corrections to the wave 
functions, we have to perform the k-space integration 
over the first Brillouin zone (BZ) in order to find the 
change in the charge density from Eq. (^. The BZ inte- 
gration is also required for calculating the incomplcte- 
basis-set corrections to the dynamical matrix. It is 
in general a full-zone integration while for the high- 
symmetry wave vectors the integrals are reduced to that 
portion of the BZ which is irreducible with respect to the 
symmetry of the perturbation vector. 

Two kinds of the integrals have to be performed in the 
linear-response calculation. The first one has the follow- 
ing form 



/i(q)=^2/k,Afe,(q) 
and the second one is given by 

/2(q) = 2/kj (l - /k±qj')^k±q/kj 

£kj ~ Ckiqj' 



(42) 



(43) 



where Akj (q) and M^^'^ are the matrix elements 
which presumably are smooth functions of wave vectors. 
In order to calculate these integrals we use the tetrahe- 
dron method in Ref. In this method, the BZ is set 
up by the reciprocal-lattice translational vectors and it 
is divided into small primitive cells exactly as in stan- 
dard fast-Fourier-transform analysis. The calculation 
becomes simpler if the q vector coincides with a mesh 
point because k± q vectors are also mesh points. In this 
way the energy bands, the expansion coefficients, and 
the structure constants have to be calculated only once 
at the mesh of the irreducible wave vectors k for the un- 
perturbed crystal. Applying symmetry operations, these 
quantities can be found for any general k. 

When applied to a semiconductor, the tetrahedron 
method is identical |-tp the special-point method of 
Monkhorst and Packed, which means that the occupa- 
tion numbers /kj in (^2|) and (^) can be regarded as the 
geometrical weights of the k-points. Both integrals (|4^ ) 
and (^) converge rapidly with respect to the number of 
k-points. The integral /2(q) reduces to the integral Ii (q) 
by performing the summation over the unoccupied bands 

For metallic systems a significantly larger number of 
k-points (Nk) is necessary when the matrix elements as 
well as the energy denominator e^j — ek±qj' are interpo- 
lated linearly within the tetrahedron. For these systems 
there are two sources of errors: the first is connected 
with the interpolation of the matrix elements and the 
second is connected with the interpolation of the Fermi 
surface. The latter can easily be circumvented in the 
linear-response calculation, since the Fermi surface can 
be determined accurately from the band structure of the 
unperturbed crystal. For the integrals /i(q) this can be 
done using the method described in Ref. ^ which is based 
on considering two, coarse and dense, meshes. In the 
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tetrahedron method the integration weight of a partic- 
ular k-point is calculated by integrating over the occu- 
pied parts of those tetrahedra that contain this point. 
The occupied part of the tetrahedron is found by linear 
interpolation between the energies at the corners of this 
tetrahedron. Suppose we introduce a much denser mesh 
that also contains the original coarse mesh. We will need 
only the energies e^j at this dense mesh, which will define 
the accurate Fermi surface (for example also by linear in- 
terpolation). Then, the occupied part of the tetrahedron 
at the coarse grid can be found by not interpolating lin- 
early the energies known at its corners but as a piece of 
the accurate Fermi surface found with help of the dense 
grid. The same is applicable to the integrals /2(q) : we 
consider the dense and the coarse grids. The band ener- 
gies are known at the dense grid. To find the integration 
weights we must find a region in the tetrahedron where 
the state |kj) is occupied and the state |k± qj) is unoc- 
cupied. This can be done using the linear interpolation 
but on the dense grid. We must also include the energy 
denominator ekj — ek±qj' ■ This is also interpolated lin- 
early but again using the dense grid. Consequently, all 
the effects from the energy bands and the Fermi surface 
are treated exactly in such scheme which allows us to 
avoid this source of errors in the integration. 

Another source of errors is connected with the linear 
interpolation of the matrix elements. We have already 
mentioned that the matrix elements are normally smooth 
functions of wave vectors and one can expect that after 
eliminating the errors connected with the approximate 
treatment of the Fermi surface, the number of k-points 
need not be too large. However, in practical calculations 
a large cancellation occurs between the two kinds of the 
integrals, ( ^^ and (43). If one uses different integration 
weights, it will lead to large numerical errors connected 
with the different convergency with respect to Nk in these 
integrals. Our task is thus to extract a large contribu- 
tion from the integral of the type /2(q) and reduce it to 
the form /i(q). This is achieved by rewriting the energy 
denominator A etj — ^k±qj' in the expression ( ^3|) as 
follows: 



1 

A 



1 



A2 + 52 



1 



A2 



(44) 



where the broadening S is usually chosen ^ O.lRy. Then, 
the sum over unoccupied bands j' in the integral con- 
taining A/(A2 + S'^) is readily performed because this 
expression remains regular when A — > 0. Consequently, 
this integral is reduced to the integral of the type Ii (q) . 
The second integral in (|4^ ) contains 1/A and must be 
treated as the integral of the type /2(q) where the orig- 
inal matrix element M^^^ is now multiplied by the 
expression in brackets in (Q). However, because the lat- 
ter rapidly goes to zero for A 3> (5, the whole integral 
remains small and it is non zero only for the band transi- 
tions j — > j' between the states near the Fermi level. In 
practical calculations of the dynamical matrix, this pro- 



cedure allows us to avoid the errors connecting with the 
large cancellations. 

We finally mention that a simple correction formula 
which significantly improves the convergency of the inte- 
grals /i(q) by taking into account the curvature of the 
matrix elements beyond the linear interpolation was de- 
rived by BlochlEllEj. Unfortunately, it is hard to derive 
such a correction for the integrals /2(q) because of the 
appearance of the energy denominator but we always use 
the Blochl correction for the integrals (^2|). 

d. Self-consistency at long wavelengths. 

The change in the charge density (^) induced by the 
displacements of nuclei screens the external perturba- 
tion (^, and the linear-response equations (|^)-(^) must, 
therefore, be solved self-consistently. Let us assume that 
we have found the response of the electrons, Sp^'^\ to the 
external perturbation SVext or the perturbation screened 
by some guessed Sp^"^^^^ (here we omit "±" for simplic- 
ity). The latter could, for instance, be the rigid shifts 
of the charge density around the displaced nuclei and in 
practical calculations the external perturbation is always 
considered as the change in the bare Coulomb potential 
(^) plus the term Vp within the MT sphere. The re- 
sponse (Jp'"-' is found along the lines described above and, 
consequently, it can be considered as some polarization 
operator H that acts on SVext , i-e. 



(45) 



If we omit the terms containing the change in the basis 
functions and forget about the completeness problem of 
the unperturbed states, the operator H is given by the 
independent-particle polarizability function tt: 



7^q(i",r') -- 



/k+qj' 



1^ fki - Ek+qj' 

xi/^k+q/(r)Vkj(r)V'k+qjv(r')V^ki(r')- (46) 

The operator tt is an integral operator while H is not nec- 
essarily one. It denotes the procedure how to construct 
the change dp from SVext ■ In particular, H contains those 
part of the operator n in which the sum over conduction 
states runs only over the number which is equal to the 
number of the basis functions, Nbasis- It also contains 
the contribution from the change in the basis functions 
according to (H), (|6l). 

After the initial response Sp^^'^ has been found, we have 
to calculate the screened perturbation (Q) . Let us call the 
Coulomb interaction, e^/jr— r'|, for vc and the exchange- 
correlation interaction in the LDA, dVxc/dp x S{y — r'), 
for Vxc- Then, the change SVeff can be written as follows: 



SVeff = SVext + [vc + Vxc)Sp, 



(47) 



and the new electronic response Sp = tiSVeff- We thus 
see that the self-consistency of the induced charge density 
means solving the Dyson equation: 
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5p = Sp'-^^ + fl{vc + v^c)Sp. 



(48) 



When q ^ 0, the integral vcSp diverges as l/<z^ which 
immediately means that searching for the solution of Eq. 
(^) by iterations, i.e. Sp = 5p^^^ + tv{vc + Vxc)5p''°^ + 
is impossible. However, it is possible when the input to 
the next, {i + l)th, iteration is prepared by mixing the 
output and input densities from the previous, (i)th, it- 
eration, i.e. Jp™*^ = XmixSpi"* + (1 - >'mix)Sp^"^, but 

the mixing parameter Xmix must be chosen to be propor- 
tional to q^. This makes the standard mixing schemes in 
the long-wavelength limit extremely time-consuming. 

This divergency problem is well known, and in the 
dielectric-matrix approach it is avoided by writing the 
solution (Esh in the form: 



Sp = e-^Sp^°\ 



(49) 



where = (1 — Tlvc — Tlvxc)^^ is an inverse dielec- 
tric matrix of the crystal. [The relation (^9|) is usually 
written for the potentials SVeff and 6Vext but in the 
present context it is more convenient to remain within 
the density language.] If for a metal live is propor- 
tional to N{ep)/q'^, where N{eF) is the density of states 
at the Fermi energy ep, then e^^ behaves as q^ when 
q — > : this is the well-known long- wavelength behav- 
ior of the metallic dielectric function. What we actually 
do when solving ( p8| ) by iterations is trying to sum up 
1 - a; + - ... = 1/(1 -f x) for a; > 1. 

In order to avoid this problem we use Thomas-Fermi- 
like screening theory. To explain the idea we assume 
that the change in charge density and the potential are 
expanded in plane waves, 



<5p(r) = EW)e 



i(q-|-G)r 



(50) 



G 



We divide the Coulomb interaction vc into a long- 
range and a short-range parts, i.e vc 



where v, 



long 
C 



lonq 

- Vc 



+ V, 



short 

c 



^■Ke^ Iq^ X ea;p[iq(r — r')]. The exchange- 
correlation in LDA is always short ranged and can be 



short 



treated together with , i.e. 

The Dyson equation can then be written as follows: 



short 



Aire 



Sp{r) = Sp^^'ir) + —^Sp{0)n^iv) + {nu;^'""'*<5p}(r), 

(51) 

where we have separated out the divergent contribution, 
Sp{0) = Sp{G = 0), and where we have called the re- 
sponse of electrons to the perturbation given by a single 
plane wave expiiqr) for nq(r). It can be written as an 
integral over the unit cell ilc- 



nq(r)= / 7rq(r,r')e^'»"-'dr'. 



(52) 



SpiO) = eY^l, (Sp^^\Q) + {n^^^^-'-SpKG = 0)) , (53) 

where we have introduced an effective dielectric constant: 
47re2 



eiong = 1 - ^nq(G = 0). 

q 



(54) 



Inserting Eq. (|5^) in to the Dyson equation ( |5l| ) we ob- 
tain: 



47re^ 



5p{t)^5p^''\y) + 

q- -T 

X {5p^°\{)) + {nw;'*''°''*(5p}(G = 0) 
xnq(r)-)-{n«;^'«"-*5p}(r). 



(55) 



The G = part of Eq. (|5l|) can be written as follows: 



where k|, = — 47re^nq(G = 0) is the Debye screening 
radius. The screened Dyson equation (|5^) is free of the 
difhculties discussed above and can be solved iteratively. 
First, one has to find the function nq(r) as the response 
of electrons to a single plane wave ea;p(zqr), and from 
that obtain Then the initial distribution 5p^^\v) is 
calculated. During the iterations one first finds the re- 
sponse to the short-range part of the perturbation, i.e. 
{IIu'*''°''*5/9}(r), and, secondly, the long-wavelength con- 
tribution is added as given by the second term in the 
right-hand side of Eq. (|5|). The output change in the 
charge density is usually mixed with the input dp to ob- 
tain an input for the new iteration. This makes the self- 
consistent cycle stable, but the mixing parameter Xmix in 
this case does not have to go to zero for q ^ and it is 
usually chosen to be 0.2 — 0.5. In practical applications 
we have found that the number of iterations required to 
solve (^5|) is about 10 while for solving the original Dyson 
equation (^ ) the number of iterations varies from 50 
to 200 depending on the length |q| of the wave vector . 
The latter is, of course, not true for those phonon modes 
where 5p{Q) = by symmetry. 

One can obviously consider the screening of not only 
the component 6V{G) - (5p(G)/|q -I- Gp with G 
but all the components within a certain sphere |q-|- G| < 
Gcutoff- This, for instance, is necessary for those zone- 
boundary wave vectors where |q| = |q -I- G|. In this case 
the function nq(r) is replaced by the functions nq-|_G(r), 
i.e. at the beginning it is necessary to calculate the re- 
sponse of the electrons to the perturbation exp[i{q+G)r]. 
The corresponding Dyson equation should be written 
again to account for the fact that e/ong is now the ma- 
trix e/ong(q+ G, q+ G'). This will reduce the number of 
iterations even more. 

Finally, we would like to point out that it should be 
possible to apply the same idea to the self-consistency 
problem in the standard band structure calculation. In 
the crystal, due to electroneutrality of the charge den- 
sity, p{G = 0) = 0. However, for those reciprocal-lattice 
vectors which are small, the components of the potential 
V{G) ~ p(G)/|G|2 might be large. This is especially 
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the case for large many-atomic unit cells. As a con- 
sequence, the mixing parameter Xmix has to be chosen 
very small. The procedure described above will require 
the calculation of the polarizability (46) with q = at 
each self-consistent iteration, i.e. the response of elec- 
trons to the plane waves ea;p[zGr] according to the ex- 
pression ( |5^ ) for all small vectors | G| < G cutoff - The 
cutoff can be chosen as the radius of the smallest first co- 
ordination sphere in G-space. The computational time 
for finding the n(5(r)-functions should presumably not 
exceed the time of one self-consistent iteration while the 
total number of iterations needed to reach the conver- 
gency is expected to decrease by approximately one order 
of magnitude, which is the case in linear-response calcu- 
lation. Note that the idea just outlined is different from 
the idea of finding the self-consistent charge transfer in 
terms of the linear-response theoryEH. For large cells we 
are screening the small G-components of the potential 
which result from some average density distribution. On 
the other hand, such details as the charge transfer be- 
tween nearest atoms is described by large G-components 
ofp(G). 



IV. APPLICATIONS. 



In recent publicationsE£[E30 we have demonstrated the 
ability of our linear-response method to compute whole 
phonon dispersions and electron-phonon interactions in 
such complicated systems as transition metals Nb and 
Mo. In the present paper we will describe application of 
the method for calculating phonon dispersions in the ma- 
terials with a few atoms per unit cell and with a relatively 
open crystalline structures. Two systems have been cho- 
sen for the applications. The first one is Si which is an 
excellent test case because of its open diamond structure. 
The second one is a transition-metal carbide NhC . This 
is a well-known classic superconductor with Tc ~ W.^K 
and its phonon dispersions show many anomalies that are 
not present in other simple-metallic and insulating sys- 
tems. Studying these anomalies as well as their influence 
on superconductivity and transport is interesting in itself 
and also represents a hard test for our method. Here we 
will only describe the calculations for the phonon disper- 
sion curves in NhC and compare the results with exper- 
iments. The calculated electron-phonon interaction and 
transport properties will be published elsewhere. 

a. Si. 



Si is a well studied elemental semiconductor from both 
experimental and theoretical sides and its jphpnon dis- 
persions have bjS^n measured Icpg time agocJ. Recent 
linear-responseQ'Q and supercelM calculations have al- 
lowed to determine its lattice dynamics for the wave vec- 
tors in the entire Brillouin zone and the results show 



a good agreement with the experiment. These calcula- 
tions were based on the linear-augmented-plane-wave 
and plane-wave pseudopotential methods. Within the 
localized-orbital representation employed in the LMTO 
method it is generally difficult to treat the materials with 
the diamond structure and, to reach close packing, inter- 
stitial spheres are usually placed into the empty sites of 
the lattice. This complicates the determination of the 
dynamical matrix. However, this problem is avoided in 
our method by the use of the Fourier transform for the 
LMTOs in the interstitial region. 

We calculate the dynamical matrix of Si as a func- 
tion of wave vector for a set of irreducible q-points in 
a (6, 6, 6)-reciprocal lattice grid (16 points per l/48th 
part of the BZ). The {I,J,K) reciprocal lattice (or 
Monkhorst-Pack) grid is defined in a usual manner: 
q.ijk = jGi + -I- -^Gs, where Gi,G2,G3 are the 
primitive translations in reciprocal space. The details 
of the calculations for every q-point are the following: 
We use 3/t — spd LMTO basis set (27 orbitals per atom) 
with the one-center expansions inside the MT-spheres 
performed up to Imax = 6. In the interstitial region, the s, 
p and d - basis functions are expanded in plane waves up 
to 15.1, 22.3, 31.7 Ry (282, 530, 868 plane waves) respec- 
tively. The induced charge densities and screened poten- 
tials are represented inside the MT-spheres by spheri- 
cal harmonics up to Imax = 6 and by plane waves with 
the 110.2 Ry energy cutoff (5208 plane waves) in the 
interstitial region. The k-integration over the BZ is per- 
formed over the (6, 6, 6) - grid (the same grid as for 
the phonon wave vectors q) by means of the improved 
tetrahedron methodS which is identical in the case of 
Si to the special-point method of Monkhorst and Pack. 
The MT-sphere radius was taken to be 2.214 a.u. and 
the Barth-Hedin-like exchange-correlation formulae af- 
ter Ref. m are employed. We use theoretically deter- 
mined lattice parameter in the calculation (the volume 
ratio V/Vexp = 0.991). 

Fig.l shows a comparison between calculated and ex- 
perimental phonon dispersion curves along the major 
high-symmetry directions. The calculated phonon den- 
sity of states is plotted at the right part in the figure. The 
theoretical frequencies are denoted by circles and the ex- 
perimental ones are denoted by triangles. The lines result 
from the interpolation between the theoretical points. 
The calculated and experimental phonon frequencies at 
the high-symmetry points F, X, and L are also listed in 
Table 1. We see that the agreement between theory and 
experiment is very good. Especially, in the optical region 
the discrepancy is about 1-1.5%, which is surprising be- 
cause the accuracy of the measured phonon modes is of 
the same order of magnitude. We also reproduce the ex- 
tended flat regions of the transverse acoustic modes indi- 
cating the accurate description of long-range interactions 
between Si -atoms as well as the correct long-wavelength 
behavior showing the good accuracy of calculated elastic 
properties of this crystal. Larger discrepancy is found for 
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the frequencies of the TA modes, where the theoretical 
branches are approximately 10% softer than the experi- 
mental ones. For instance, the calculated frequency of the 
mode is 4.00 THz while Wexp(^TA) = 4.49 ± 0.06 



THa23. The same kind of discrepancy has also been re- 
cently reported in Refs. | 16, The agreement is slightly 



improved when we recalculate the dynamical matrix at 
the X - point using the experimental lattice constant. 
We have found that the frequency of the Xta ~ mode 
is increased from 4.00 to 4.27 THz. This shows that the 
mode has a large negative Gruneizen parameter and it 
is thus very sensitive to the unit-cell volume used in the 
calculation. Because of the large LMTO basis sets, large 
I max and plane-wave energy cutoffs this discrepancy is 
hard to relate to the internal parameters in the calcula- 
tion and it is most likely connected with the local density 
approximation. 

b. NbC. 

The lattice dynamical properties of transition-metal 
carbides and, especially, NbC have attracted much atten- 
tion in the past because of the existence of pronounced 
anomalies in its acoustic branches and their influence to 
superconductivity . While some model calculations of 
the phonon dispersions exist in the literature and var- 
ious meclianisms explaining these anomalies have been 
proposedc3, no ab initio investigation of the lattice dy- 
namics for NbC have so far been performed. Here we ap- 
ply the linear-response approach to the phonon spectrum 
of NbC in order to check the accuracy of our method. 

The dynamical matrix of NbC is calculated at the 29 
irreducible q - points of a (8, 8, 8) reciprocal-lattice grid. 
The self-consistent calculations performed for every wave 
vector involve the following parameters: 3k— spd LMTO 
basis per Nb atom (27 orbitals) and 3k — sp LMTO ba- 
sis per carbon atom (12 orbitals). The one-center ex- 
pansions inside the MT-spheres are performed up to 
Imax = 6. In the interstitial region the basis functions 
are expanded in plane waves up to 13.4, 19.6, 26.9 Ry 
(136, 228, 338 plane waves) for, respectively, s,p and d - 
orbitals of Nb, and up to 24.1, 35.8 Ry (306, 536 plane 
waves) for s,p -orbitals of C. The changes in the charge 
densities and the potentials are represented inside the 
MT-spheres by spherical harmonics up to Imax = 6 and 
by plane waves with an 121 Ry energy cutoff (3382 plane 
waves) in the interstitial region. The k-space integra- 
tion for the matrix elements is performed over a (8, 8, 8) 
- grid (the same grid as for the phonon wave iiectors q) by 
means of the improved tetrahedron methodLJ. However, 
the integration weights for the k-points of this grid have 
been found to take into account the effects arising from 
the Fermi surface and the energy bands precisely. This 
is done with help of a (32,32,32) grid (897 k - points 
per 1/48 BZ) as we explained in Section III(c) of this 
paper. The MT-sphere radius of Nb is taken to be 2.411 
a.u. and the radius of the carbon sphere is 1.786 a.u. 



The Barth-Hedin-like exchange-correlation formulae af- 
ter Ref. 39 are employed. As in the case of Si, we also 
use the theoretically determined lattice parameter in this 
calculation (the volume ratio V/Vexp = 0.982). 

The results of our calculations are presented in Fig. 
2, where we compare theoretically determined phonon 
dispersions (circles) with those measured by inelastic- 
neutron-scattering techniquecil (triangles). The calcu- 
lated phonon density of states is plotted at the right part 
of the figure. The lines are simply the result of inter- 
polation between the theoretical points. Since the q— 
grid (8, 8, 8) considered here is still too coarse to resolve 
the anomaly of the longitudinal acoustic branch near the 
wave vector (0.6, 0, 0) 27r/a , we have performed a sepa- 
rate calculation for the q -point (0.625, 0, 0) which fits to 
the (16, 16, 16) - grid in k-space. We see that the agree- 
ment between theory and experiment is good. Most of 
the calculated frequencies agree within a few percent with 
those measured despite of the fact that we have used only 
29 k-points for the BZ-integration. (We list for compar- 
ison our calculated and experimental phonon frequencies 
at the high-symmetry points F, X, and L in Table 2.) 
The theory reproduces the major anomalies presented in 
the acoustic branches: the well-known anomaly near the 
wave vector (0.6, 0, 0) 2n/a which is also present and well 
describedE3 within our linear-response method in pure 
Nb crystal; the anomaly of the longitudinal mode near 
the wave vector (0.5,0.5,0) 27r/a as well as large soften- 
ing of both TA and LA modes near the L -point. More- 
over, we also predict an anomalous behavior of the lowest 
transverse acoustic mode along the (^^0) direction near 
the wave vector (0.5, 0.5, 0) 27r/a .Here the frequencies 
are not known experimentally. The anomaly found by 
us is, however, less pronounced compared to-thc results 
of double-shell model calculations of WebercS while we 
have certainly not too many points along this direction 
to judge about its exact dispersion. 



V. CONCLUSION. 

In conclusion, we have described in detail an all- 
electron linear-response approach based on density func- 
tional theory and the LMTO technique. The method 
is developed to calculate lattice dynamical properties of 
crystalline solids and is uniquely applicable for the sys- 
tems with broad and narrow energy bands. For test pur- 
poses, we have applied the method to compute phonon 
dispersions for Si and NbC which have open structures 
and two atoms per unit cell. The results of our appli- 
cations are in a good agreement with the experiment. 
We have thus shown that accurate calculations of lat- 
tice dynamics are now possible even for such compli- 
cated systems as jtjjansition-metal compounds. In the 
forthcoming paperE^ we give a description of our method 
for calculating electron-phonon interactions and apply 
the method to compute lattice-dynamical, superconduct- 
ing and transport properties for a large number of ele- 
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mental pCnetals (a brief report of this work has appeared 
alreadyll3) . In another pubhcationC3 we describe an ap- 
phcation of the method for computing electron-phonon- 
coupHng strength in Ca — Sr — Cu — O high- Tc super- 
conductor. 
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APPENDIX: Fourier transform of pseudoLM- 
TOs. 

Consider a Hankel function i^'„i(r) = K K,i{r)i''Yim{Y) 
of energy which is singular at the origin. The three- 
dimensional Fourier transform of this function JCkL (k) is 
known to behave as fc'~^ for large k. The task is to substi- 
tute the divergent part of K^i{r) inside some sphere s by 
a smooth regular but otherwise arbitrary function. This 
function is chosen so that the Fourier transform is con- 
vergent fast. In the full-potential LMTO method of Ref. 
p8| , the augmenting function is the linear combination of 
the Bessel function J^l and its energy derivative J^l 
matched together with its first-order radial derivative 
with the Hankel function at the sphere boundary. The 
Fourier transform becomes convergent as k~'^. One can 

(«) 

obviously include higher-order energy derivatives J kL in 
order to have a smooth matching at the sphere up to the 
order n. This was done in connection with the problem 
of solving the Poisson equation in Ref. |3^. The Fourier 
transform here converges as fc~(3+n) ^■\^^ prefactor 
increases as {21 4- 2ri -I- 3)!! and this prohibits the use of 
large values of n. A similar procedure has been also used 
in the LMTO method of Ref. In the present work we 
will use a different approach based on the Ewald method. 
Instead of substituting the divergent part only for r < s 
we consider the solution of the equation: 

(56) 

The function on the right-hand side of the Helmholtz 
equation is a decaying Gaussian. The parameter a/ is a 
normalization constant: a; — ■\/277r(2?7^)'+''/^s^'+^/(2Z — 
1)!!. The most important parameter is rj. It is chosen 
such that the Gaussian is approximately zero when r > s 
and rj must depend on / as well as the sphere radius s. 
The solution K^lIy) is thus the Hankel function for large 
r, it is a regular function for small r and it is smooth to- 
gether with its radial derivatives at any r. The function 
Ki^i^r) can be calculated in terms of the following error- 
function-like contour integral: 



(57) 



When 77 — cxD this integral is known as the Hankel in- 
tegral. The most important result is that the Fourier 
transform of K^i (r) decays exponentially. It is given by: 



9 J+i /-oo 



(58) 



Restoring the original notations, the pseudoLMTOs 
XkJ?x i^) ^^'^ tti^ Bloch waves of wave vector k as defined 
in Eq. (p7|). The Fourier coefficients XkrlO^ + G) are 
given by: 



47r s 



i+i 



|k + G| 



X.RLi^ + G) - ^^(2Z-1)!! |k + G|2-K2 
g(K^-|k+GP)/4,Ly^(k + G)e-'(''+°)^, (59) 

where flc is the volume of the unit cell and where we have 
subscripted rj with the indexes Rl and s with R. 

In practical calculations the parameter rjm can be 
chosen from the ratio between the Hankel function 
at the sphere and the solution of Eq. (p6|), i.e. 
K ki{s r) / K t^i{s b) = 1 + 5. The error \5\ is usually taken 
not larger than 0.03 which leads to the number of plane 
waves per atom needed for the convergency in (p7| ) vary- 
ing from 150 to 250 when I ~ 2. For the s,p— orbitals this 
number is smaller by a factor of 2 — 3. 
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TABLE I. Comparison between calculated and experimental phonon frequencies at the high-symmetry points F, X, and L 
for Si [THz]. 





Flto 


Xta 


Xlao 


Xto 


Lta 


Lla 


Lto 


Llo 


theory 


15.56 


4.00 


12.27 


13.90 


3.09 


11.20 


14.78 


12.38 


exp". 


15.53 


4.49 


12.32 


13.90 


3.43 


11.35 


14.68 


12.60 



"Reference y 



TABLE II. Comparison between calculated and experimental phonon frequencies at the high-symmetry points F, X, and L 
for NbC [THz]. 



Flto 


Xta 


Xla 


Xto 


Xlo 


Lta 


Lla 


Lto 


Llo 


theory 17.05 


6.37 


7.51 


17.64 


18.65 


4.26 


6.02 


18.82 


21.60 


exp". 16.70 


6.35 


7.30 


17.20 


17.80 


4.00 


6.00 




19.20 



"Reference [41| 
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Fig. 1. Calculated phonon dispersion for Si (circles) along the high-symmetry directions in comparison with the 
experimental (triangles). The lines are the result of interpolation between theoretical points. Also shown is the 
calculated phonon density of states (DOS). 

Fig. 2. Calculated phonon dispersion for NbC (circles) along the high-symmetry directions in comparison with 
the experimentcil (triangles). The lines are the result of interpolation between theoretical points. Also shown is the 
calculated phonon density of states (DOS). 
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